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Abstract — This paper illustrates the application of recent re- 
search in region-of-attraction analysis for nonlinear hybrid limit 
cycles. Three example systems are analyzed in detail: the van 
der Pol oscillator, the "rimless wheel", and the "compass gait", 
the latter two being simplified models of underactuated walking 
robots. The method used involves decomposition of the dynamics 
about the target cycle into tangential and transverse components, 
and a search for a Lyapunov function in the transverse dynamics 
using sum-of-squares analysis (semidefinite programming). Each 
example illuminates different aspects of the procedure, including 
optimization of transversal surfaces, the handling of impact maps, 
optimization of the Lyapunov function, and orbitally-stabilizing 
control design. 

I. Introduction 

The purpose of this paper is to illustrate a new technique 
for estimation of regions of attraction for nonlinear hybrid 
limit cycles proposed in (U. Three example systems have been 
chosen which illuminate different aspects of the method. 

A major motivation for the work in this paper is control of 
underactuated "dynamic walking" robots (|2|). These robots 
can exhibit efficient, naturalistic, and highly dynamic gaits. 
However, control design and stability analysis for such robots 
is a challenging task since their dynamics are intrinsically 
hybrid and highly nonlinear. 

Estimates of regions of attraction can be useful for many 
problems, including analysis of different candidate control 
laws, generation of tree-based feedback motion-planning con- 
trollers (0), or planning transitions among a library of pre- 
stabilized walking gaits ((H). 

The method involves choosing a decomposition of the 
dynamics into tangential and transversal components, and then 
searching for a Lyapunov function in the transversal compo- 
nents that verifies a "tube" about the limit cycle in which 
orbital stability is guaranteed. The verification is performed 
using sum-of-squares (SoS) programming. 

The first example is the van der Pol oscillator, chosen 
because it is very well-known and well studied, and thus 
provides a good test of the method. With this example, we 
illustrate the importance of selecting the transversal decompo- 
sition intelligently. Indeed, it is shown that the commonly-used 
orthogonal transversal surfaces are often a bad choice. 

The second is the "rimless wheel", a one-degree-of-freedom 
hybrid mechanical system, which serves as a simple model of a 
walking robot. The dynamics of the sytem are simple enough 
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that much can be said about the rimless wheel analytically. 
Despite being simple, it exhibits hybrid (switching) behaviour 
representing the foot fall of a walking robot. 

The third example is the "compass-gait" walker, a more 
complex model of a walking robot. For this system, one cannot 
derive analytical regions of attraction, so a computational 
approach is essential. This system requires transversal surface 
optimization, proper handling of impacts, and also orbitally- 
stabilizing control design based on a transverse linearization. 

A. Background 

The most well-known tool for analysing limit cycles is the 
Poincare map: orbital stability is characterized by stability of 
an associated "first-return map", describing the repeated passes 
of the system through a single transversal hypersurface. Often 
a linearization of the first-return map is computed numerically, 
and its eigenvalues can be used to verify local orbital stability. 
Since the system's evolution is only analyzed on a single 
surface, regions of stability in the full state-space are difficult 
to evaluate via the Poincare map. 

A related technique known variably as "transverse coordi- 
nates" or "moving Poincare sections" also has a long history 
and was certainly known to exist by Poincare, however has not 
been much used in applications until recently due to difficulty 
in the relevant computations (see Q). With this technique, a 
new coordinate system is defined on a family of transversal 
hypersurfaces which move about the orbit under study. In most 
cases, it is also used to study local stability, however it can 
be adapted to characterize regions of stability in the full state 
space. 

The method we propose is to construct the transverse 
dynamics in regions of the orbit, and then utilize the well- 
known sum-of-squares (SoS) relaxation of polynomial positiv- 
ity which is amenable to efficient computation via semidefinite 
programming (see, e.g., (6), Q, 0). The sum-of-squares 
relaxation has been previously used to characterize regions 
of stability of equilibria of nonlinear systems (see, e.g., l9l . 
fTOlL EB) and as a tool for constructive control design in 0. 

There is comparatively little work on computing regions 
of stability of limit cycles. The proposed method has aspects 
in common with the surface Lyapunov functions proposed in 
C21, however that method was restricted to piecewise linear 
systems. The technique of cell-to-cell mapping, proposed by 
fT3lL improves the efficiency of exhaustive grid-based methods 
of regional analysis and has been used in analysis of simple 



walking robots ( lfT4l ). however the computational cost is still 
exponential in the dimension of the system. 

II. Problem Statement 

We consider the following class of hybrid systems with 
planar switching surfaces: 



x = f(x, u), x ^ S 
c + = AO), x e S~. 



(1) 

(2) 



where x G M n , and u G M m . Suppose /(•) and A(-) are 
smooth and A : S~ —> where 



S = {x : c'_x = 

£ + — {x : c' + x = <i + }, 



(3) 
(4) 



c_,c + G M n , and G R. Suppose #*(•) is a non- 

trivial T-periodic solution that undergoes N impacts at times 
{ti, t2, tjsf} + fcT for integer fc. We will assume that 
the impacts are not "grazing", i.e. c f _f(x*(ti)) ^ and 
c' + /0r*(^)) ^ for alii. 

It is well-known that periodic solutions of autonomous 
systems cannot be asymptotically stable, since perturbations 
in phase are persistent. The more appropriate notion is orbital 
stability (see, e.g., 0, (H), EH)- 

The analysis objective is to efficiently compute a region of 
state space D C W 1 from which orbital stability to #*(•) is 
guaranteed. 

III. Transverse Dynamics and Regions of Orbital 
Stability 

In this section we briefly describe the process of verifying 
regions of orbital stability. Full details of each step are given 
in Q. 

The process we propose for finding regions of orbital 
stability is based on the construction of a smooth local change 
of coordinates x — » (x±,t). At each point t G [0, T] we define 
a hyperplane S(t), with 5(0) = S(T), which is transversal to 
the solution #*(•), i.e. x*(t) S(t). 

The transversal surfaces are defined by: 

S(T) = {y£R n :z(Ty(y-x*(T))=0} 

where z : [0, T] — >> R n is a vector function which will be 
optimized. 

Given a point x nearby £*(•), the scalar r G [0, T) 
represents which of these transversal surfaces S(r) the current 
state x inhabits; the vector x± G R n_1 is the "transversal" 
state representing the location of x within the hyperplane S(r), 
with x± = implying that x = x*(r). This is visualised in 
Figure [T] 

The process for computing regions of orbital stability is as 
follows: 

1) Select a family of transversal surfaces S(r), and the 
associated transformation x — > (x±,r) such that at im- 
pact times the transversal surfaces line up with switching 
surfaces. 

2) Compute the nonlinear dynamics in this new coordinate 
system as well as a periodic linear system representing 





V{x±,t) < 1 



Fig. 1. Top: transversal surfaces S(r) around the target orbit x*, with 
a particular solution x(t) converging to the orbit x*. Bottom: a Lyapunov 
function defined on a transversal surface. 



the dynamics of close to the orbit: the transverse 
linearization. 

3) Construct a candidate quadratic Lyapunov function as- 
sociated with the transverse linearization via standard 
techniques from linear control theory. 

4) Using this result as an initial seed, iteratively solve 
a sequence of sum-of-squares programs to compute 
maximal regions in which a Lyapunov function can be 
found verifying both well-posedness of the change of 
coordinates and orbital stability for the true nonlinear 
dynamics. 

For each r G [0,T], U(r) G R(n-i)xn is a smoot h matr i x 
function of r projecting x ^ x_\_. 

The dynamics in the new coordinates x±,r are given by 
the continuous dynamics 



x ± 



= f *pv 

n(r)/0r*(r))f, 

z(TYf(x*(T)+U(T)'x ± ) 



U(r)'x ± + U(r)f(x*(r) + U(t)'x ± ) 

(5) 
(6) 



z(t)7(x*(t))-^'iI(t)'s ± ' 
for t^U, and impact updates 

x+ =n(r+)[A,( a; *(rr)+n(rr)'x:) -x*(t+)], (7) 

when t = t{. The change of coordinates and the above 
dynamics are well-defined in a region around the target orbit 

£*(•). 
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A. Regions of Orbital Stability 

Suppose there exists a Lyapunov function V such that 
V(x±,t) > 0,x± ^ 0,V(0,r) = for all r e [0,T] 
for which the following conditions hold on the level set 

{x:V{x ± ,t)<1}: 

d :V(x ± ,T) < -5\x ± \ 2 , 



dt 



z(rYf(x*(T))-^-'u(TYx ± > 0, 



(8) 
(9) 



for some S > 0. The first constraint verifies stability via 
the decreasing Lyapunov function; the second verifies well- 
posedness of the change of variables. If hybrid dynamics 
are considered, at the switching times one must verify the 
condition: 

y(n(r+)[A i (x*(rr)+n(rr)V) 

-x*(t+)],t+)-V(x ± ,t-) < 0, (10) 

If these conditions hold then the set {x : V(x±,r) < 1} is 
an inner (conservative) estimate of the region of attraction to 
the limit cycle. Furthermore, for polynomial dynamics these 
can be relaxed to a SoS program. In some cases, one must 
also check that the impact surface is not reached before it is 
expected (see Q). 

Typically we sample a sufficiently fine finite sequence of 
re [0, T] and verify the conditions on x± at each sample. In 
the optimization procedure, one must search for both V(x±,r) 
and Lagrange multipliers verifying the regional conditions. 
This problem is non-convex, however when fixing V(x±,r) 
and searching over Lagrange multipliers it is a semidefinite 
program, and when fixing Lagrange multipliers and searching 
over V(x±,t) it is a semidefinite program. Thus if one has 
a reasonable initial guess for V(x±,r) an iterative procedure 
can be applied. Whilst this is not guaranteed to find a global 
optimum, in practice the authors have found it works very 
well. 

B. Transverse Linearization 

In the construction of initial candidate Lyapunov functions 
and feedback controllers, we will make use of the linearization 
of the transverse dynamics: 



x± 



A(t)x ± (t) 
A d x± : t~- 



■B(t)u(t), t + U 



(11) 
(12) 



representing the parts of §5§ and ([7]) linear in x±. Note that 
this can be given analytically for any system in the class. 

IV. The van der Pol Oscillator 

The first example we consider is the van der Pol oscillator, 
defined by the following differential equation: 

y- v(i-y 2 )y + y = o. 

where p is a constant, which we take equal to 1. It is well- 
known that this system has a single unstable equilibrium at 
the origin, and a periodic cycle which is the limit from every 




Fig. 2. Verified regions of orbital stability using constant rescalings of linear 
Lyapunov function and orthogonal transversal surfaces. 



other point in the plane. Thus, it makes a good simple example 
on which to test the proposed method. The above differential 
equation can obviously be rewritten in terms of a state x = 
[x\ X2\' := [y y)' in the form 



x = f(x), f(x) 



(1 



x 2 

x\)x 2 - Xi 



To compute the periodic solution x* we simply simulated the 
system forward from an initial condition away from the origin 
until it converged. Let T be its period. 

Given a vector z(r), for a planar system the projection 
II(r) onto surfaces orthogonal to z(r) is simple: II (r) = 
[—Z2(r) Zi(t)]. We start with transversal surfaces orthogonal 
to the system motion, i.e. z(r) = f(x*(T))/\f(x*(r))\. 

A natural candidate for a Lyapunov function is the solu- 
tion of the Lyapunov differential equation for the transverse 
linearization: 



P(t) + A(t)'P(t) + P(t)A(t) + Q = 0. 



(13) 



For any Q > a unique periodic solution P(t) = P(t + T) > 
exists, and suggests a Lyapunov function of x±P(r)x±, as 
was suggested in fT5ll . 

We can then search for the maximal level set x±P(r)x_\_ < 
p in which stability can be verified. In the framework of 
Section III-A| this corresponds to verifying ^ and ^ for 
Lyapunov fuction of the form V(x) = (l/p)x±P(r)x±, 
which can be performed via a simple bisection search over 
p. The results are shown in Figure [2] 

The regions are strongly limited by points at which the 
change of variables x — » (x±,r) becomes ill-defined. In the 
figure, these points resemble the hub of a bicycle wheel. Math- 
ematically, they correspond to points at which the denominator 
of ([6]) goes to zero. This is clearly a consequence of the choice 
of transversal surfaces and motivates exploring other possible 
choices. 

The van der Pol oscillator is a special case, since we know 
in advance that an orbital stability test must fail at the origin, 
the unstable equilibrium. With this in mind we construct radial 
transversal surfaces centred at the origin. Since z(r) and II (r) 
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Fig. 3. Verified regions of orbital stability using constant rescalings of linear 
Lyapunov function and radial transversal surfaces. 




Fig. 4. Verified regions of orbital stability using time-varying cr(r) and 
radial transversal surfaces. 

are changed, we recompute the solution of ( [T3] ) and again 
search for the maximal level set of this Lyapunov function. 

The results are plotted in Figure [3] We no longer have the 
problem with singularities, however the regions are still quite 
thin. The reason is that whilst the Lyapunov ODE is guaranteed 
to locally verify orbital stability (see (HI), it may not be a very 
good choice for verifying regional stability. 

The result is much improved if we allow time-varying 
adjustment to the Lyapunov function via a scaling function 
cr(r) > 0, i.e. 

V(x ± ,t)=*(t)x' ± P(t)x ± (14) 

where a(-) : [0,T] R, cr(0) = <j(T), and P(r) is as above. 

The results are plotted in Figure [4] This figure is computed 
with a(r) a Bezier polynomial of order 20. It is observed that 
increasing the order of a(r) grows the region essentially to 
the origin. 

Note that since we are searching over a class of candidate 
Lyapunov functions which are symmetric with respect to the 




Fig. 5. The function cr(r) used to verify the region in Fig El 




Fig. 6. Verified regions of orbital stability using time-varying cr(r) and 
locally-optimized transversal surfaces. 

orbit, the computed region is the best that can be achieved. The 
time-varying rescaling that achieved this is shown in Figure 

13 

To obtain these results we made use of qualitative knowl- 
edge of the true region of stability in order to choose the 
transversal surfaces, which will not be possible for more 
complicated systems. In [1] a procedure for choosing z(r) 
was suggested which depends only on local information about 
f(x*,u*). The idea is to optimize the function z(r) so that 
the distance to the closest point of singularity is maximized. 
The results of applying this procedure are shown in Figure [6] 

This choice of transversal surfaces does not do quite as well 
the radial surfaces, which is not surprising. However it does 
substantially better than the orthogonal surfaces in Fig. [2] It 
seems there is still plenty of room for improvement in the 
optimization of z(r), perhaps in an iterative procedure based 
on results of a prior region-of- stability computations. 

V. Rimless Wheel 

The rimless wheel is a simple planar model of walking 
and consists of a central mass with several 'spokes' extending 
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Fig. 7. Rimless wheel system 

radially outward. (See Figure [7]). At any given moment one of 
the spokes is pinned at the ground, and the system follows the 
dynamics of a simple pendulum, f{0,6) = [0, sin(0)]'. When 
another spoke contacts the ground, the system undergoes an 
inelastic collision governed by + = cos(2a)0, and the new 
spoke becomes the pinned one. 

On a sufficiently inclined slope the system has a stable 
limit cycle, for which the energy lost in collision is perfectly 
compensated by the change in potential energy. The rimless 
wheel has been analyzed in the literature and the basin of at- 
traction has been computed exactly (see ifTTl ). It is interesting 
to see how close our region of attraction estimation method 
can approximate the actual basin for this system. 

Figure [8] shows the phase portrait of the rimless wheel, 
with arrows indicating the direction of the dynamics. The right 
edge of the graph represents the collision surface that maps 
to the left edge of the graph (or vice-versa, depending on the 
direction of dynamics). Because the impact depends only on 
the value of the angle, 9, and not on 9, the collision surfaces 
are vertical. The thick green lines are the homoclinic orbits of 
the simple pendulum. The thick black line shows the stable 
limit cycle, and the shading shows a subset of its true region 
of attraction. 

In this case, it is natural to select vertical transversal 
surfaces, since there are no singularities in the change of 
variables, and the transversal surfaces are aligned with switch- 
ing surfaces. The surfaces can be parametrized by 9 as 
t = t(9), and the nominal trajectory is simply 9* (9) = 9 and 
9* (0) = ^2E - 2cos(0), where E is the total energy of the 
system. Then the transversal coordinate is the vertical position 
with respect to the nominal trajectory: x± — 9 — 6* (6). 

With this choice of transversal surfaces, the dynamics in the 
new coordinate system are straightforward to compute: 



x ± 



dt 
sin(0) 



0*(0) = sin(0) 
_sm(0) - = _ 

6*(0) 



_ d9*_ 

~ ~d9 
sin(0) 



(9)9 



"(0) 



■x±. 



(15) 



Interestingly, the transversal dynamics are linear for any given 
9. Note that the transverse dynamics are also easily derived 
from equations ([5J and ([6J by setting z(r) = [1,0]' and 
II(r) = [0,1]. This gives f = an d the equation ( p3] ), 

as expected. 




0.1 0.2 0.3 
(radians) 

Fig. 8. Phase portrait of the rimless wheel system. 
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Fig. 9. Regions of Attraction for the rimless wheel limit cycle. Light shaded 
region is the true RoA. Dark shaded region is the verified RoA. 



As with the van der Pol oscillator, to find an initial candidate 
Lyapunov function we computed the unique periodic solution 
of the jump Lyapunov differential equation: 

- P(t) = A(t)'P(t) + P(t)A(t) + Q, t + U 
P(r-) = A d (T z yP(T+)A d (r z ) + Q z , t = U 

with Q, Qi > and search over scalar rescalings: V(x) = 
(1/ P )x> ± P(t)x ± . 

Figure [9] shows the results. The discrete set of transverse 
surfaces can be seen as thin black vertical lines around the 
orbit. The computed basin of attraction (dark shading) is 
within the true basin (light shading), but doesn't fill it fully. 
This is not surprising, since we searched over a very restrictive 
set of Lyapunov functions. 



6 




0.1 0.2 0.3 

(radians) 



Fig. 10. Regions of Attraction for the rimless wheel limit cycle. Light shaded 
region is the true RoA. Dark shaded region is the verified RoA. 




Fig. 11. Compass gait system 



We then searched over 10 th -order rescaling polynomials 
cr(r), as in ( [T4| ), to maximize area of the computed region. 
Figure [10| shows the results and, as expected, the basin of 
attraction is larger. 

Note that the verified basin of attraction includes regions 
where 9 < and, by the choice of transverse surfaces, f < 0. 
The verification procedure still holds when the system state 
moves backwards through the transverse surfaces. 

VI. Compass-Gait Walker 

The compass-gait walker is a two-degree-of-freedom (four- 
state) nonlinear hybrid system shown in Figure [TT] Similarly to 
the rimless wheel, one of the legs (the "stance leg") is always 
pinned at the ground, and the other (the "swing leg") swings 
until it hits the ground, at which point the swing leg becomes 
the stance leg and vice versa. A motor at the hip generates 
torque, r, between the legs. Letting q = [6 sw ,6 st ] f , u = r, 
and / = a + b the continuous dynamics can be expressed in 
the standard manipulator form as 

H(q)q + C(q,q)q + G(q)=Bu, 



where 
H = 

C = 

G 



mb 2 

-mlbcos(6 st - 6 SW ) 


mlbsm(6 st - 6 SW )6 SW 
mbg sm(6 sw ) 
■(rrihl + rna + ml)gsm(0 st ) 



J SW ) 

,2 



- mlbcos(9 st 
(rrih + m)l 2 + ma 

mlbs'm(6 st - sw )0 st 


" 1 



-1 



At the moments of impact, the coordinates undergo a relabling: 
®fw — ®st and 6f t = 6 SW . The impact dynamics for velocities 
are derived assuming a perfect inelastic collision: 

QU + = QaQ~ 

where 

-mab + (rrihl 2 + 2mal) cos(2a) 
—mab 



-mab 




mb(b - 



-Zcos(2a)) ml (I — b cos(2a)) + ma 2 
mb 2 —mbl cos(2a) 



m h f 



and a = Usw ~ Ust . The impact takes place when SW + st + 
27 = 0, where 7 is the slope of the ground. Note that although 
the switching dynamics are highly nonlinear, the switching 
surfaces are planar, matching the assumptions of Q. 

For some combination of parameters and initial conditions 
the system exhibits a limit cycle behavior of walking downhill 
passively (with zero torque). For this reason, the model is often 
used for studying bipedal walking (see, e.g. fl8l . fT9lL l20lL 
EH , l22l and many others). Prior to this work, studies of 
the small basin of attraction of this system were limited to 
exhaustive simulation. 

We orbitally stabilize the limit-cycle using LQR control in 
the transverse coordinates, and analyze the region of attraction 
for the closed loop system. We begin by computing the T- 
periodic stable limit cycle solution numerically from suitable 
initial conditions. We next optimize for a set of transversal 
surfaces at N = 40 sample points along the trajectory. These 
surfaces are constrained to align with planar switching surfaces 
capturing the "swing leg" losing then regaining contact with 
the ground. We construct the transverse coordinates, and the 
resulting transverse dynamics. 

A set of transversal surfaces were chosen via the optimiza- 
tion procedure suggested in [1]. An LQR controller was then 
designed for the transverse linearization by solving for the 
periodic positive-definite solution of the jump-Ricatti equation: 

-P = A'P + SA- PBR- X B'P + Q, t^U 
P(rr) = A d (Ti)'P(Tt)A d (Ti) + Qi, t = U 

with R, Q, Qi > 0. The feedback control is given by: 

u(t,x ± ) = -R- 1 B(t)'P(t)x ± . 

Note that there is no nominal control command, since we are 
stabilizing a passive walking cycle. 

To enable the use of SoS optimization, we approximate 
the non-polynomial dynamics of the walker via a third order 
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Taylor expansion. This expansion is performed around each of 
the N sample points. 

The Riccati equation provides us with an initial candidate 
Lyapunov function: 

Vq(x±,t) = a x f ± P(r)x±. 

Taking a suitably small cr , we find a valid region of attraction. 
To optimize the region of attraction, we rescale Vo by a 
polynomial p : [0, T] \-> (0, oo): 



V(x ± ,t) 



1 



p( r , 



■V (x±,t). 



We maximize the integral of p(t) over the interval as a 
surrogate for the volume of the region of attraction. For 
the compass gait system this parameterization proved better 
numerically conditioned than directly scaling the Lyapunov 
function. The following observations allow our constraints to 
remain linear in the parameters of p(t). Note that: 



d m( \ 1 



p(t)^V (x ± ,t) - ^tVo(x±,t) 



To enforce -^V(x±^r) < — S\x±\ 2 , we can exploit the fact 
that S can be an arbitrary positive constant, and instead require: 



p(t)^Vo(x±,t) 



d_ 

dr 



p(r)fV (x ± ,r) < S \x±\ 



As p{t) is bounded above on [0,T] an appropriate 5 exists. 
Finally, {x \ V(x±,t) < 1} is identically {x | Vq(x±,t) < 
p(r)}. Figure 12 presents overlapped plots of the region of 
attraction projected into the and (62,62) planes. 

These regions indicate that the controller can stabilize much 
larger variations in the swing-leg than in the stance-leg, which 
matches intuition on the compass gait walker. 



A. Accuracy of Taylor Expansion 

Since we are approximating the dynamics by a Taylor 
expansion, it is important to check whether we can trust 
the verification. To examine this approximation, we sample 
points on the boundary of the certified region of attraction. For 
each point, we compute ^V(xj_^r) with the true dynamics 
and Taylor expansion approximation. At each time step we 
sampled 10,000 points from the boundary. In Figure [13] we 
plot the true value of 



j^V against the approximation error. 



Note that the vertical scale of the plot is approximately 10 -5 
vs. an almost unit scaling along the horizontal, and essentially 
all the points are to the left of zero. 

Of the 400,000 samples, 104 had f t V > 0, so the verifica- 
tion using Taylor series was correct at 99.97% of samples. The 
largest value of j^V was 0.0026. Note that isolated samples 
of positive ^ V do not necessarily imply instability from that 
point, it is possible (indeed, likely) that the Lyapunov function 
was not appropriate for proving stability at that point. 



Compass Gait Projected Region of Attraction 




.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 



Fig. 12. Regions of Attraction for the compass-gait walker limit cycle. We 
present two of six possible projections of the RoA. The upper curve plots the 
trajectory and RoA projected into the (61,61) plane, that of the "swing leg". 
The lower curve is the equivalent for the "stance" leg coordinates (#2, 02 )• 



dV/dt Error due to Taylor Expansion 




True dV/dt 



Fig. 13. The above compares the rate of change of the Lyapunov function for 
the true and Taylor expanded dynamics at the boundary of the approximate 
region-of-attraction. 99.7% of samples were verified correctly. 



B. Computation Times 

In computing the regions of attraction, we alternate between 
optimizing over V(x±,r) - a "V-step", and optimizing over 
Lagrange multipliers - an "L-step". It took four V-L iterations 
until covergence. The computations were performed on a 
2.66 GHz intel Core i7 processor with 8GB of RAM, and 



timings for each iteration are shown in Table VI-B The total 
computation time was around 17.7 minutes. 

Note, however, that due to the sampling in r, each L-step 
is made up of 41 independent optimizations (40 continuous 
samples and one impact map) each of which took between 
1.6 and 2.6 seconds. Since these computations are trivially 
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Iteration 


1 


2 


3 


4 


L-step time (s) 


75.0 


90.2 


80.8 


78.7 


V-step time (s) 


194 


183 


201 


170 



TABLE I 

Computation times for the compass-gait walker. 



parallelizable, substantial speed increases could be achieved 
on multiprocessor machines. 

VII. Summary and Future Work 

The purpose of this paper has been to demonstrate the ap- 
plication of a new method for estimating regions of attraction 
for nonlinear hybrid limit cycles. Three simple examples were 
chosen that elucidated important aspects of the technique, 
including selection of transversal surfaces, the handling of 
impacts, optimization of the Lyapunov function, and control 
design. 

This work can be extended in a number of directions. It 
will serve as an essential component in recent algorithms for 
feedback control and motion planning in [3]. Implementation 
on more complex models and experimental validation will be 
an important test. 

In this paper we searched for regions via r-varying rescal- 
ings of quadratic Lyapunov functions from linear theory. This 
was done primarily to simplify the optimization procedure, 
however in principle it is possible to search over any poly- 
nomial Lyapunov functions, which would presumably give 
substantially bigger (and possibly asymmetric) regions for 
many systems. 

There are many connections between sum-of- squares ver- 
ification and robust control and stability theory via integral 
quadratic constraints ( 1231 , |24|) which will be investigated. 
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